    function [stforc,yforc]=zb_forecast_loop(GG,RR,ZZ,CC,sT,Emat,nforc) 
    % ====================================================================
    % function [stforc,yforc]=zb_forecast_loop(GG,RR,ZZ,CC,sT,Emat,nforc) 
    % sT   [ns 1] vector of last state 
    % Emat [nforc nx] matrix of shocks if empty RR is not used 
    % nforc is number of forecasts
    % AJ 1/18/2011 
    % =====================================================================
    nst=length(sT(:));
     
    if ~isempty(Emat)
        dim=size(Emat);
        if dim(1)~=nforc
            error('EMAT Matrix of shocks must be [NFORC NX]')
        end
        flag_useShocks=1;
    else
        flag_useShocks=0;
    end
    
    % Define matrices     
    stforc=zeros(nst,nforc);
    yforc =zeros(size(ZZ,1),nforc);

    % Fill in first state 
    switch flag_useShocks 
        case 0
            stforc(:,1)=GG*sT(:);
        case 1 
            stforc(:,1)=GG*sT(:) + RR*( Emat(1,:)' ) ;
    end 
    yforc( :,1)=ZZ*(CC+stforc(:,1) );
    for ii=2:nforc;
        switch flag_useShocks
            case 0
                stforc(:,ii)=GG*stforc(:,ii-1);
            case 1
                stforc(:,ii)=GG*stforc(:,ii-1) + RR*( Emat(ii,:)' ) ;
        end 
        yforc(:, ii)=ZZ*(CC+stforc(:,ii) );
    end
    stforc=stforc';
    yforc =yforc' ;